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Abstract 

Long-lived electronic coherences in various photosynthetic complexes at cryogenic and room tempera- 
ture have generated vigorous efforts both in theory and experiment to understand their origins and explore 
their potential role to biological function. The ultrafast signals resulting from the experiments that show 
evidence for these coherences result from many contributions to the molecular polarization. Quantum pro- 
cess tomography (QPT) is a technique whose goal is that of obtaining the time-evolution of all the density 
matrix elements based on a designed set of experiments with different preparation and measurements. The 
QPT procedure was conceived in the context of quantum information processing to characterize and under- 
stand general quantum evolution of controllable quantum systems, for example while carrying out quantum 
computational tasks. We introduce our QPT method for ultrafast experiments, and as an illustrative exam- 
ple, apply it to a simulation of a two-chromophore subsystem of the Fenna-Matthews-Olson photosynthetic 
complex, which was recently shown to have long-lived quantum coherences. Our Fenna-Matthews-Olson 
model is constructed using an atomistic approach to extract relevant parameters for the simulation of pho- 
tosynthetic complexes that consists of a quantum mechanics/molecular mechanics approach combined with 
molecular dynamics and the use of state-of-the-art quantum master equations. We provide a set of methods 
that allow for quantifying the role of quantum coherence, dephasing, relaxation and other elementary pro- 
cesses in energy transfer efficiency in photosynthetic complexes, based on the information obtained from 
the atomistic simulations, or, using QPT, directly from the experiment. The ultimate goal of the combina- 
tion of this diverse set of methodologies is to provide a reliable way of quantifying the role of long-lived 
quantum coherences and obtain atomistic insight of their causes. 

* Contributed equally to this work. These three authors are ordered in the sequence that their main scientific contri- 
bution is emphasized in the text. 



aspuru @ chemistry.harvard.edu 



1 



I. INTRODUCTION 



The initial step in photosynthesis is highly efficient excitonic transport of the energy captured 
from photons to a reaction center [IJ. In most plants and photosynthetic organisms this process 
occurs in light-harvesting complexes which are interacting chlorophyll molecules embedded in a 
solvent and a protein environment Several recent experiments show that excitonic coherence 
can persist for several hundreds of femtoseconds even at physiological temperature [l3]-[6l. These 
experiments suggest the hypothesis that quantum coherence is biologically relevant for photosyn- 
thesis. The results have motivated a sizeable amount of recent theoretical work regarding the 
reasons for the long-lived coherences and their role to the function. 

The focus of many studies is on the theoretical models employed. In this context, it is essential 
to be as realistic an possible and employ the least amount of approximations. Most of the currently- 
employed methods involve a master equation for the reduced excitonic density operator where the 
vibrational degrees of freedom (phonons) of the protein and solvent are averaged out. Amongst 
these simple methods are the Haken-Strobl model and Redfield theory as employed in Refs. [|71[8]| 
and [HI respectively. To interpolate between the usual weak and strong exciton-phonon coupling 
limits, Ishizaki and Fleming developed a hierarchical equation of motion (HEOM) theory which 
takes into account non-equilibrium molecular reorganization effects [fTOll . Jang et al. perform 
a second order time-convolutionless expansion after a small polaron transformation to include 
strong coupling effects [fTTTl .Another set of studies focuses on the role of quantum coherence and 
the phonon environmentin terms of transport efficiency or entanglement. It was shown that the 
transport efficiency is enhanced by the interaction or interplay of the quantum evolution with the 
phononic environment lEHll [121 • Entanglement between molecules is found to persist for long 
times [[T3UT51. 

The ongoing effort can be summarized with two equally important questions: What are the 
microscopic reasons for the persistence of quantum coherence and what is the relevance of the 
quantum effect to the biological functionality of the organism under study? In this work, we sum- 
marize the recent efforts from our group to approach the problem from several angles. Firstly, 
we investigate the role of coherences in the exciton transfer process of the Fenna-Matthews-Olson 
(FMO) complex. We quantify the amount and the contribution of coherence to the efficient energy 
transfer process. Secondly, we present our quantum mechanics/molecular mechanics (QM/MM) 
approach to obtain information about the system at the atomistic level, such as detailed bath dy- 
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namics and spectral densities. Finally, we propose a spectroscopic tool that allows for obtaining 
directly the information of the quantum process via our recent theoretical proposal for the quantum 
process tomography technique to the ultrafast regime. 

II. THE ROLE OF QUANTUM COHERENCE 

In this section, we discuss the question about the relevance of quantum effects to the bio- 
logical function. A negative answer to this question would mean that a particular effect, while 
being quantum, is not leading to any improvement in the functionality of a biological system, and 
therefore would be a byproduct of the spatial and temporal scales and physical properties of the 
problem. For example, in energy transfer (FT) quantum coherence could arise from the closely 
packed arrangement of the chromophores in a protein scaffold but it could, in principle, represent 
a byproduct of that arrangement and not a relevant feature. Another example, it may be true that 
the human eye can detect a single photon, but it is not clear if this quantum effect is relevant to the 
biological function, which usually operates at much larger photon fluxes. If, on the other hand, the 
above yes-no question of the relevance is answered positively for a particular effect in a biological 
system, it would present a major step towards establishing the relevance or importance of a quan- 
tum biological phenomenon. A natural follow-up questions is: How important quantitatively is a 
particular quantum effect? 

Both of these questions should preferably be studied by experimental means. An experiment 
would have to be designed in a way that tests for the biological relevance of quantum coherence. 
Possible experiments could involve quantum measurements on mutated samples. In the FMO 
complex that acts as a molecular FT wire the efficiency of the transport event is most likely a good 
quantifier for biological function. One would need a way to experimentally quantify this efficiency 
and extract the relevance of quantum coherence to the efficiency. This can be hard in practice. Yet, 
as we will discuss in this work, quantum process tomography is able to obtain detailed information 
about quantum coherence and the phonon environment and might thus lead to progress in this area. 

In the case when experimental access to an observable that involves the biological relevance 
is hard or impossible, a theoretical treatment can provide insight. It is illustrative to analyze a 
model of the particular biological process in terms of a quantifier for the success of the process. 
An example is the aforementioned efficiency of energy transport. In bird vision, the quantum yield 
of a chemical reaction is a relevant measure [16,1 . Once a detailed model and a success criterion 
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is established, one needs to quantify the contribution of quantum coherence to the success crite- 
rion. For this step, one can proceed in two distinct pathways. The first pathway is a comparison 
to a classical reference point; the success criterion is computed for the actual system/model and a 
classical reference model that does not include quantum correlations. The difference of these two 
values is attributed to quantum mechanics and can be considered the quantum mechanical contri- 
bution to the success of the process. For example, the energy transfer dynamics of a sophisticated 
quantum mechanical model such as [[TOl could be compared to a semi-classical Forster treatment 
that leads to a hopping description. In general, this comparison strategy has the drawback that one 
has to invoke a classical, and in some cases very artificial, model. 

Our work has been mainly concentrated on a second theoretical pathway in answering the 
relevance question, which overcomes this issue. It is based on just the quantum mechanical model 
and the success quantifier. No other, for example classical, model is invoked. The actual model will 
contain dynamical processes that are quantum coherent and others that are incoherent. The non- 
trivial task is to deconstruct how the various processes contribute to the performance criterion. This 
can be done by decomposing the performance criterion into a sum of contributions, each associated 
with a particular process. The terms in this sum related to quantum mechanical processes will then 
give a theoretical answer to the overall relevance of the particular process and will quantify this 
relevance. This line of thought was developed and discussed in Ref. [fTTl for energy transfer in the 
FMO complex and provided insight into both questions "Is a quantum effect relevant?" and "If yes, 
how much?", at least from a theoretical standpoint within the approximations of the model under 
consideration. In this section, we extend this idea to include the effect of the initial conditions and 
compare the results to a total integrated coherence, or concurrence, measure. We utilize secular 
Redfield theory and the hierarchy equation of motion approach. 

The Hamiltonian describing a single exciton is given by: 

He = J^(em + X)\m){m\ + Jmn i\m){n\ + \n){m\) . (1) 

m m<n 

where the site energies e^, and couplings Jmn are usually obtained from detailed quantum chem- 
istry studies and/or fitting of experimental spectra. The reorganization energy A, which we assume 
to be the same for each site, is the energy difference of the non-equilibrium phonon state after 
Franck-Condon excitation and the excited-state equilibrium phonon state. The set of states \m) 
is called the site basis and the set of states \a) with He\a) = Ea\a) is called the exciton basis. 
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We now briefly introduce the secular Redfield master equation in the weak exciton-phonon (or 
system-bath) coupling limit and the non-perturbative hierarchy equation of motion approach. In 
both approaches, the dynamics of a single exciton is governed by a master equation, which is 
schematically given by: 

d 

—pit) = Mpit) = {Mn + Aldecoherence + A^trap + -Mloss) pit). (2) 

The master equation consists of the superoperator Ai, which is divided into several components. 
First, coherent evolution with the excitonic Hamiltonian is described by the superoperator 
Ain = —i[He, ■]. In addition, decoherence due to the interaction with the phonon bath is incor- 
porated by decoherence- decoherence dcpcuds ou the spcctral dcusity, which modcls the coupling 
strengths of the phonon modes to the system. Finally, one has the processes for trapping to a 
reaction center A^trap and exciton loss A^ioss due to spontaneous emission. Associated with these 
processes are the trapping rate n and the loss rate F. Details about the trapping and exciton loss 
processes can be found in ifTTlfTSll . 

The secular Redfield theory is valid in the regime of weak system-bath coupling. The superop- 
erator decoherence IS of Liudblad foHTi with Liudblad operators for relaxation in the exciton basis 
and for dephasing of excitonic superpositions. The relaxation rates depend on the spectral density 
evaluated at the particular excitonic transition frequencies, satisfy detailed balance, and depend on 
temperature through the Bose-Einstein distribution. The dephasing rates are linear in temperature. 
We use the same Ohmic spectral density as in [fTOl , i.e. J{u) = 2\'^uo/Tt{uj'^ + 7^), where I/7 is 
the bath correlation time. For I/7 = 50 fs, this spectral density shows only modest differences to 
the spectral density used in ifTTll . Further details about the Lindblad model can be found in [fTTl . 

The hierarchy equation of motion approach IfTOl consistently interpolates between weak and 
strong system bath coupling. The assumption that the fluctuations are Gaussian makes the second- 
order cumulant expansion exact. The resulting equation of motion can be expressed as an infinite 
hierarchy of system, i.e. and connected auxiliary density operators {cTj}, arranged in tiers. 
For numerical simulation, "far-away" tiers in the hierarchy are truncated in a sensible manner. 
The hierarchy equation of motion can also be written as in Eq. ([2]) when we make the replacement 
p{t) — )■ (p(t), (Ti, 0-2, ■ ■ ■ ) and use the hierarchical structure discussed in flOl for the decoherence 
superoperator Al decoherence- For simulations of the Fenna-Matthews-Olson complex, we use the 
scaled hierarchy approach developed in [.19,1 . It was shown recently that four tiers of auxiliary 
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density operators are enough for accurate room temperature simulations [|20l . which enables the 
rapid computation of efficiency and total coherences. The trapping and exciton loss processes are 
naturally extended to the auxiliary systems. 



In our previous work [fTTll . we developed a method to quantify the role of quantum coherence to 
the transfer efficiency. The energy transfer efficiency (ETE) is given by the integrated probability 
of leaving the system from the sites that are connected to the trap instead to being lost to the 
environment. That is, 77 = dtTi{Aitra.pPit)}- It was shown that the ETE can be partitioned 
into r] = r]ii + r/decoherencej whcrc the efficiency due to the coherent dynamics with the excitonic 
Hamiltonian is given by: 



The ETE contribution r^decoherence involves A^dccoherence, i-e. r^dccoherence = Tr{ A^trap ( A^trap + 

A^ioss)^^-^dccoherence-^^p(O)}. lu tMs work, wc cxtcud our ETE Contribution method to quan- 
tify the role of the initial state to the ETE. We obtain a separation of the coherent contribution, 
Vh = ^init + Vdyn, whcrc the efficiency r^init can be ascribed to the initial state. The T^dyn is defined 
by T^dyn = Vn — ^init and can be interpreted as dynamical part of the coherence contribution arising 
during the time evolution. For the computation of r^init, we note that one can always express the en- 
semble described by the system density matrix as p{t) = Pinit(i)|^imt(i))(V^init(^)|+X]fcPfc(^)Pfc('^)- 
Here, Pinit(i^) is the probability of the quantum system being in the (Hamiltonian time-evolved) 
initial state |V^init(^))5 where J5init(0) = 1. The Pkit) are the probabilities of being in some other 
ensemble state pk{t), where Pinit(i) + YlikPk{t) = 1- The probability Pinit(^) is modified by 
the interaction with the environment and readily computed for Markovian Lindblad dynamics by 
considering the damped no-jump evolution due to the decoherence superoperator Aldecoherence 
ifTSl I2TI l22l . Therefore, we can compute the efficiency pertaining to the initial state by r^init = 
dtTT{Mtra.pPmit{t)\i^mitit)){'^imt{t)\}. Togcthcr with Equatiou ([3]), this obtains the desired 
separation ?7h = Vmit + Vdyn- 

Additionally, we employ another measure for the role of coherence by straightforwardly inte- 
grating over time all the coherence elements of the density matrix. That is: 



r/H = Tr{A^ trap (A^ trap + -Mioss)-'A^ H^^ " V(0) } • 



(3) 




00 



dt I {m\ p (t) \n) 



(4) 
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We normalize with respect to the case of coherent evolution at A = 0.0/cm, i.e. C{\) = 
C(A)/C(0). Based on the discussion in [fT3l . the quantity C can be considered as the (nor- 
malized) integrated entanglement (concurrence) that is present before the exciton is trapped in the 
reaction center or lost to the environment. We note that the total coherence measure C is similar in 
spirit to a measure of the first kind discussed above. This is because the normalization essentially 
performs a comparison of the actual model at a certain A with an artificial model at A = 0. (For 
the numerical evalutation, the integral in Eq. Q is computed until Tr{p(t)} < 10^'^.) 




Reorganization energy [ 1 /cm] Reorganization [ 1/cm] Reorganization [ 1 /cm] 



FIG. 1 . (Left panel) Efficiency rj (solid black) and contributions of initial state r/init (dash-dotted gray) and 
coherent evolution r^init + 7?dyn (dashed red) for a dimer that is based on the strongly coupled sites 1 and 2 
of the Fenna-]VIatthews-OIson complex using the secular Redfield model. The initial state is at site 1 and 
the target is site 2. At a physiological value of around A = 35/cm, one finds r/init = 0.0 and ?7dyn = 0.43. 
(Center panel) Efficiency and integrated coherence C for the dimer with the secular Redfield approach. At 
A = 35/cm there is C = 0.37. (Right panel) Same quantities as in the center panel for the dimer using 
the hierarchy equation of motion approach with 15 tiers of auxiliary systems. At A = 35/cm, one finds 
C = 0.44. The parameters are 1/k = 1 ps, l/F = 1 ns, and I/7 = 50 fs for all panels. 

In Fig. [1} we present the two measures of coherence for a dimer system. For the dimer, we 
take the sites 1 and 2 of the FMO complex with ei = 0/cm, 62 = 120/cm, and J = — 87.7/cm, 
see [|23l . and room temperature. This system will also be the focus of the following sections on 
the atomistic detail simulations and quantum process tomography. Here, for studying the role of 
quantum coherence, we assume that the task is defined by the exciton initially being at the lower 
energy site 1 and the target site being site 2. In the left panel of Fig. [T] we show the efficiency 
1], the contribution r/n from Eq. (|3]), and r/init for the secular Redfield model. In the present small 
system, environment-assisted transport is relatively unimportant, with the efficiency as a function 
of the reorganziation energy being close to unity everywhere. The underlying contributions show 
a transition from a regime dominated by coherent evolution to a regime dominated by incoherent 
Lindblad jumps. At A = 35/cm, we find r/init = 0% and ?/h = 43%. In Fig. [T] (center panel), 
we find that the total coherence measure C for the dimer is around 0.37 for A = 35/cm. In Fig. 
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FIG. 2. (Left panel) Efficiency rj (solid black) and contributions of initial state r/init (dash-dotted gray) 
and coherent evolution rji^it + ^dyn (dashed red) for the Fenna-Matthews-Olson complex using the secular 
Redfield model. The initial state is a classical mixture of site 1 and 6 and the target site for trapping is site 3. 
The actual system has a reorganization energy of around A = 35/cm, where ?7init = 0.0 and r/dyn = 0.17. 
(Center panel) Efficiency for initial site 1 (solid black) and initial site 6 (dashed black) and integrated 
coherence C for initial site 1 (dashed red) and initial site 6 (dash-dotted green) for the Fenna-Matthews- 
Olson complex with the secular Redfield approach. At A = 35/cm there is C = 0.0151 (inital site 1) and 
C = 0.0017 (initial site 6). (Right panel) Same quantities as in the center panel for the FMO complex using 
the scaled hierarchy equation of motion approach with four tiers of auxiliary systems. At A = 35/cm, one 
finds C = 0.020 (inital site 1) and C = 0.0022 (initial site 6). The parameters are 1/k = 1 ps, l/T = 1 ns, 
and 1/7 = 50 f s for all plots. 

[T] (right panel), the total coherence is plotted for the dimer in the hierarchy equation of motion 
approach. We use 15 tiers of auxiliary systems. At A = 35/cm, we find C = 0.44; because of the 
sluggish, non-equilibrium bath there is more coherence than in the secular Redfield model. 

In Fig. [2] (left panel), we present the coherent, decoherent, and initial state contribution to 
the ETE for the Fenna-Matthews-Olson complex as a function of the reorganization energy for 
the secular Redfield model at room temperature. We use the Hamiltonian given in [12311 and the 
contribution measures given in Equation ([3]) and by r^init. The initial state is a classical mixture 
of site 1 and 6. For small reorganization energy, the efficiency is around t] = 60% and for larger 
reorganization energies we observe environment-assisted quantum transport (ENAQT) Q, with 
the efficiency rising up to almost i] = 100% for the physiological value of A = 35/cm. The 
contributions measures T^dyn and T^init reveal the underlying dynamics. The quantum dynamical 
contribution r/dyn is around 17% at A = 35/cm In our model, this part is due to an interplay of 
the Hamiltonian dynamics and the trapping/loss dynamics, which both have their preferred basis 
being the site basis. The main part of the efficiency at A = 35/cm is due to incoherent Lindblad 
jumps, having a value of r/decoherence = 83%. The initial state contribution is relevant only at small 
values of the reorganization energy. 

' In Ref. ifm . we found the value 7/h = 10% for a different Hamiltonian and a different spectral density. 
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In Fig. [2] (center and right panel), we compare the efficiency and the coherence measure C for 
the secular Redfield and the hierarchy equation of motion approach [fTOl for the Fenna-Matthews- 
Olson complex. The initial state is either localized at site 1 or at site 6. Four tiers of auxiliary 
systems were used in the computation, which already lead to a good agreement with [TOl for the 
dynamics at A = 35cm^^, I/7 = 50 fs, and room temperature. In Fig. [2](right panel), ENAQT is 
observed with increasing reorganization energy also in the hierarchy approach, with the efficiency 
rising up to almost i] = 100% at A = 35/cm. In Fig. |2] (center and right panel), it is observed 
that the normalized total coherences of the density matrix decrease with increasing reorganization 
energy. For the secular Redfield case, we obtain C(\ = 35cm^^) = 0.0151 for the initial site 
1 and C(\ = 35cm^^) = 0.0017 for the initial site 6. For the hierarchy case, we obtain more 
coherence, i.e. ^(A = 35cm-^) = 0.020 for the initial site 1 and ^(A = 35cm-^) = 0.0022 for 
the initial site 6. In both models, coherence is more important for the rugged energy landscape of 
the pathway from site 1 than for the funnel-type energy landscape of the pathway from site 6. 

Master equation approaches, such as the ones discussed in this section suffer from various 
drawbacks. Redfield theory is only applicable in the limit of weak system bath coupling and does 
not take into account non-equilibrium molecular reorganization effects. The hierarchy equation of 
motion approach assumes Gaussian fluctuations and Ohmic Drude-Lorentz spectral densities. The 
detailed atomistic structure of the protein and the chlorophylls is not taken into account in these 
approaches. The results thus provide a general indication of the behavior of the actual system but 
not a conclusive and detailed theoretical proof. In the next section, we will present a first step 
toward such a detailed study with our combined molecular dynamics/quantum chemistry method. 
The atomistic structure is included and realistic spectral densities can be obtained. We also present 
a straightforward method to simulate exciton dynamics beyond master equations. We thus address 
the second question of the microscopic origins of the long-lived quantum coherence. 



III. MOLECULAR DYNAMICS SIMULATIONS 

Among many other biologically functional components, protein complexes are essential com- 
ponents of the photosynthetic system. Proteins remain as one of the main topics of biophysical 
research due to their diverse and unidentified structure-function relationship. Many biological 
units are highly optimized and efficient, so that even a point mutation of a single amino acid in 
conserved region often results in the loss of the functionahty. [|24] - l26l Have the photosynthetic 
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system adopted quantum mechanics to improve its efficiency in its course of evolution? To answer 
this question, careful characterization of the protein environment to the atomistic detail is neces- 
sary to identify the microscopic origin of the long-lived quantum coherence. As explained in the 
previous section, the contribution of the quantum coherence to the energy transfer efficiency in 
biological systems have been successfully carried out, yet a more detailed description of the bath 
in atomic detail would be desirable to investigate the structure-function relationship of the pro- 
tein complex and to test validity of the assumptions used in popular models of the photosynthetic 
system. 

The site energy of a chromophore is a complex function of the configuration of the chro- 
mophore molecule, and the relative orientation of the molecule to that of the embedding protein 
and that of other chromophore molecules. Factors affecting site energies have intractably large 
degrees of freedom, so it is reasonable to treat those degrees of freedom as the bath of an open 
quantum system. The state of the system is assumed to be restricted to the single exciton man- 
ifold. To construct a system-bath relationship with atomistic detail of the bath, we start from 
the total Hamiltonian operator, and decomposed the operator in such a way that the system-bath 
Hamiltonian is not assumed to be any specific functional form: 

Htotal ^^em{Rch,Rprot)\m){m\ +y^^{Jmn{Rch,Rprot)\m) {n\ + C.C.} 

m m,n 

~l~ '^ch ~l~ Tprot ~l~ ^chip' 1 Rchi Rprot) ~l~ ^prot^Rchj Rprot)- 

e^n represents the site energy of mth site, J„in is the coupling constant between mth and nth sites. 
cr denotes the excitonic state of chromophores, Rch corresponds to the nuclear coordinates of 
chromophore molecules, and Rprot are the nuclear coordinates of the remaining protein and en- 
closing water molecules. T and V are the corresponding kinetic and potential energy operators for 
the chromophores and proteins respectively under Bom-Oppenheimer approximation. The poten- 
tial energy term for chromophores depends on the exciton state of the systen, because dynamics 
of a molecule will be governed by different Born-Oppenheimer surface when its excitonic state 
changes. However, as a first approximation, we assumed that the change of Born-Oppenheimer 
surfaces does not affect the bath dynamics significantly. With this assumption, we can ignore the 
dependence of the excitonic state in the Vch term and the system-bath Hamiltonian only contains 
the one way influence from the bath to the system: 
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Htotal ~ y^^em{Rch,Rprot)\m){m\ +y^^{Jmn{Rch,Rprot)\m) (n| + C.C.} 

m m,n 

+ ^m{Rch, Rprot) \^) ("^| + T^h + Tpj-ot + K/il-Rc/i, Rprot) + ^rot (-Rc/i, -R; 

m 

= ^e„|m)(m| + ^ { Jmn|"^)(n| + c.c.} 



+ ^ {em(ilcfc, -Rprot) - em} 1"^) ("^1 + 



^ ^ \^Jmn{R'ch^ Rprot) Jmn\ -|- C.C. 



~l~ -^c/i ~l~ Tpyot ~\~ Vchi^R'ch: R'prot) ~\~ Vproti^R'chi R'prot] 



Hb 

(6) 

Based on this decomposition of the total Hamiltonian, we set up a model of the FMO com- 
plex in atomistic detail with the AMBER force field [|27l |28|| and approximate the propagation 
of the entire complex by classical mechanics. Molecular dynamics simulations were conducted 
at 77K and 300K with an isothermal-isobaric (NPT) ensemble. The parameters for the system 
and the system-bath Hamiltonian were calculated using quantum chemistry methods along the tra- 
jectory from the molecular dynamics simulation, was calculated using the Q-Chem quantum 
chemistry package. Il29l The electronic excitations were modeled using the time-dependent den- 
sity functional theory using the Tamm-Dancoff approximation. The density functional employed 
was BLYP and the basis set employed was 3-2 IG*. External charges from the force field were 
included in the calculation as the electrostatic external potential. The coupling terms, Jm„, were 
obtained from the Hamiltonian presented in [|23l and considered to be constant in time. e„i was 
chosen as time averaged site energy for the mth site to minimize the magnitude of the system- 
bath Hamiltonian. In this work, only site 1 and site 2 were considered for the exciton dynamics. 
However, the methodology can be applied for the exciton dynamic of all seven chromophores. 

To obtain a closed-form equation for the reduced density matrix, we applied mean-field ap- 
proximation ll30l : because no feedback from the system to the bath was assumed, the state of the 
bath is not affected by the state of the system. Therefore, the total density matrix, W{t), can be 
factorized into the reduced density matrix p{t), and B{t) which is defined only in the Hilbert space 
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FIG. 3. (Left panel) Time evolution of the exciton population at the site 1 (pu) based on the strongly coupled 
site 1 and 2 of the FMO complex at 77K and 300K. The initial pure state p = |1)(1| was propagated using 
Monte Carlo integration of unitary evolutions, where the time-dependent site energies are obtained from a 
combined molecular dyanmics/quantum chemistry approach. The asymptotic distribution does not follow 
a Boltzmann distribution because relaxation of the system to the bath is not considered. (Right panel) The 
concurrence between site 1 and 2 at 77K and 300K. Quantum coherence lives longer at a lower temperature. 

of the bath. With additional assumption that the bath is in thermal equilibrium, we can obtain the 
closed equation for the reduced density matrix. 



= 'l [Hs,p{t)] - '-Tr {[HsB,Wm 

^ \Hs,p{t)] - [TT{HsBB{t)},p{t)] (7) 

Thermal equilibrium of the bath was ensured by the thermostat of the molecular dynamics sim- 
ulation. Thus, the reduced density matrix was obtained by Monte Carlo integration of 4000 in- 
dependent instances of unitary quantum evolution with respect to the thermally equilibrated bath. 
Each instance was propagated by integrating the Schrodinger equation with the simple exponential 
integrator. 

Figure [3] shows the change of the population of the site 1, pn, and the concurrence between 
site 1 and 2. The population is evenly distributed between the two sites because relaxation was not 
considered. The concurrence, 2|pi2|, is an indicator of pairwise entanglement for the system. [fT3ll 
Note that the coherence builds up during the first ^ 100 fs , and then decreases subsequently due 
to the decoherence from the bath. 

Figure |4] shows the spectral density of the first chromophore. Although the spectral density 
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FIG. 4. (Left panel) Spectral density from the autocorrelation function of the site 1 of the FMO com- 
plex from the molecular dynamics simulation at 77K and 300K. While the spectral density reflects the 
characteristic vibrational modes of the protein and the chromophore molecule, high-frequency modes are 
overpopulated due to the limitation of the Newtonian mechanics. (Right panel) Absorption spectrum of site 
1 and 2 at 77K and 300K. 



of the bath from molecular dynamics simulation shows characteristic frequencies related to the 
actual protein environment and the bacteriochloropyll molecule, high-frequency modes are over- 
populated due to the limitation of the classical mechanics. There are efforts to incorporate quantum 
effects into the classical molecular dynamics simulation in a slightly different context, [13114331 and 
we are investigating the possibilities of applying these corrections. 

Another simplification employed was the omission of the feedback from exciton states. When 
the exciton state of a bacteriochlorophyll is changed, the Bom-Oppenheimer surface which gov- 
erns the dynamics of the chromophore molecule should be also changed. The current scheme only 
propagates the protein complex on the electronically ground-state surface. Incorporating the feed- 
back could lead to the different characteristics of the protein bath. There exist several schemes for 
mixed quantum-classical dynamics [|34] - l36l which potentially resolve the problem at the additional 
computational cost of simultaneously propagating excitons and protein bath. 

Calculations are underway to carry out the full seven-site simulation of the FMO complex at 
different temperatures to compare with experimental temperature-dependent results 

In the following final section, we will describe our quantum process tomography scheme, which 
is a spectroscopic technique associated with a computational procedure for direct extraction of the 
parameters related to the quantum evolution of the system, in terms of quantum process maps. 
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IV. QUANTUM PROCESS TOMOGRAPHY 



So far, we have delved into several theoretical models to characterize quantum coherence in 
the entire FMO complex and in a dimer subsystem of it. Experimentally, however, a clear char- 
acterization of this coherence is still elusive. Signatures of long lived quantum superpositions 
between excitonic states in multichromophoric systems are potentially monitored through four 
wave-mixing techniques |i3l|37l|38l. However, a transparent description of the evolving quantum 
state of the probed system is not necessarily obtained from a single realization of such experiments. 
In these, a series of three weak incoming ultrashort pulses sent from a noncoUinear setup induce 
a macroscopic third order polarization in the sample. The latter manifests in a time dependent 
spatial grating through which a fourth pulse diffracts. From an operational standpoint, this last 
pulse selects the spatial Fourier component of the polarization which corresponds to its wavevec- 
tor (heterodyne detection), hence earning the name of four wave-mixing for this technique (FWM) 
ll38ll . Extracting specific Fourier components of the induced polarization allows for the selection 
of a particular set of processes in the density matrix of the probed system, as each wavevector is 
associated with a carrier frequency of the pulse. These processes can be intuitively understood 
by keeping track of the dual Feynman diagrams that account for the perturbations that the pulses 
induce on the bra or ket sides of the density matrix of the probed system. Whereas the analysis of 
these experiments is naturally carried out in the density matrix formalism, an important question 
is whether the density matrix itself can be imaged via these experiments, a problem known as 
quantum state tomography (QST) [|39l . If this were possible, quantum process tomography (QPT) 
could also be carried out, therefore providing a complete characterization of excited state dynamics 
Boll . In a previous study, we showed that a series of two-color heterodyned rephasing photon-echo 
(PE) experiments repeated in different polarization configurations yields the necessary informa- 
tion to carry out QST and QPT of the single-exciton manifold of a coupled heterodimer [41]. In 
the present article, we adapt our previous theory to extract this information from two-dimensional 
spectra, similar to those employed in current experiments. 

We begin by reviewing some basic aspects of QPT. Under very general assumptions, the evo- 
lution of an open quantum system can be described by a linear transformation [l42|: 




(8) 



cd 
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where PabiT) is the element ah of the reduced density matrix p of the system at time T. Equation 
^ is remarkable in that x{T) is independent of the initial state. Knowledge of x{T) implies a 
complete characterization of the dynamics of the reduced system and, in fact, QPT can be opera- 
tionally defined as the procedure to obtain x{T). Conceptually, it is straightforward to recognize 
that, due to linearity, x{T) can be inverted by preparing a complete set of inputs, evolving them 
for time T, and detecting the outputs along a complete basis. In the context of nonlinear optical 
spectroscopy, this is exactly the strategy we shall follow, with a few caveats due to experimental 
constraints. 

To place the discussion in context, we shall be again concerned with the subsystem composed 
of the excitonic dimer between sites 1 and 2 of the FMO complex. For simplicity, we ignore 
the rest of the sites in this theoretical study. We only need to be concerned with four eigenstates 
of this model system: The ground state \g), the delocalized single-excitons \a) and and the 
biexciton |/), which in the photosynthetic system can be safely assumed to be the direct sum of 
the single-excitons without significant interactions between them. Therefore, the biexciton energy 
level is just ujf = Ua + oj/s- We label the delocalized excitons so that |a) is the higher energy 
eigenstate compared to \(3). Denoting the transition energies between the i-th and the j-th states 
by Uij = LOi — tUj, it follows that tUag = tOf/3 and tUfSg = ujfa [l43l . The excitonic system is not 
isolated, and in fact, it interacts with a phonon and photon bath which induces relaxation and 
dephasing processes in it. 

The experimental technique we consider is photon-echo (PE) spectroscopy, which is a par- 
ticular subset of FWM techniques where the wavevector of the fourth pulse corresponds to the 
phase-matching condition kpE = — fci + ^2 + ks, with fcj being the wavevector corresponding to 
the i-th pulse. Here, the labeling of the pulses corresponds to the order in which the fields interact 
with the sample. Typically, the ultrashort pulses employed to study these excitonic systems possess 
an optical carrier frequency, therefore allowing transitions which are resonant with the frequency 
components iw^g and ±uJag- In PE experiments, the first pulse centered at ti creates an optical 
coherence beating at a frequency Uga or w^^. At ^2 = ^1 + t, the second pulse creates a coherence 
or a population in the single exciton manifold. At ^3 = ^2 + T, the third pulse generates another 
optical coherence, but this time, beating at the frequencies opposite to the ones in the first interval, 
that is, at frequencies coag or co^g, causing a rephasing echo of the signal. The heterodyne detection 
of the nonlinear polarization signal Ppe{t, T, t) occurs at time ^4 = ^3 + 1. Borrowing from NMR 
jargon, the intervals (ti, ^2), (^2, ^3), and (^3, ^4) are traditionally refered to as coherence, waiting, 
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and echo times, and their durations are r, T, and t, respectively. This nomenclature should not be 
taken literally. For example, in most case, coherences do not only evolve in the coherence time, 
but in the waiting and echo times. Similarly, the waiting time is often referred to as population 
time, which hosts dynamics of both populations and coherences. For a historical perspective on 
this vocabulary, we refer the reader to any comprehensive NMR treatise such as [|44|. 

The experiment is systematically repeated for many durations for each interval. In order to 
'watch' single-exciton dynamics, it is convenient to isolate the changes on the signal due to the 
waiting time T. This exercise is accomplished by performing a double Fourier transform of the 
signal along the r and t axes, which yields a 2D spectra that evolves in T [|45l - l47ll : 



In order to map a PE experiment to a QPT, we identify the coherence interval as the prepa- 
ration step and the echo interval as the detection step. This assumption implies that the optical 
coherence intervals have well characterized dynamics. This hypothesis is reasonable due to a sep- 
aration of timescales where optical coherences will presumably decay exponentially due to pure 
dephasing and not due to intricate phonon-induced processes. Therefore, the 2D spectrum consists 
of four Lorentzian peaks centered about {ur, oJt) = {ujag, ^ag), {^ag, ^pg), {^Pg, ^ag), {^i3g, ^Pg)- 
In this discussion, we shall ignore inhomogeneous broadening, noting that it can always be ac- 
counted for as a convolution of the signal with the distrubtion of inhomogeneity. The width of 
these Lorentzians can be directly related to the dephasing rates of the optical coherences. Loosely 
speaking, a particular value on the cUj- axis of the spectrum indicates a specific type of state prepa- 
ration, whereas the Ut axis is related to a particular detection. More precisely, a peak in the 
2D spectrum displays the correlations between the frequency beats from the coherence and echo 
intervals. A crucial realization is that the amplitude of these peaks can be written as a linear com- 
bination of elements of the time evolving excitonic density matrix stemming from different initial 
states, that is, of elements of x{T) itself [1411 : 




(9) 
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S{ujo,g,T,Uag) = • ei)(/u„p • 62) 

x{C(^3[(Mag • e3)(Ma3 " e4)(X99aa(T) - 1 - Xaaaa(T)) 
+ (M//3 • e3)(M//3 • e4)X/3^aa(r)] 

+Cf3[(M/a • e^){^l'f0 ■ 64) - (M/3s • e3)(Ma5 " e4))Xa/3aa(r)]} 
^{CZ^[{^^ag ' e3){tlag ' e4){XggfSa{T) (T)) 

• e3)(M//3 • e4)xmaiT)] 

+C^s[{ifJ'fa ■ e3)(M//3 • 64) - • 63) (Mas " e4))Xa/3/3a(r)]}, (10) 



S(uJag,T,Uj0g) = -C^,C^^{llag ' ei)iHag ' 63) 

(XOLOiO. iT)] 

+CZM^^fP ■ e3)(M/a • 64) - (Mas • e3)(M/3s • e4))x^aaa(7')]} 

-dCf,(/^a<,-ei)(M/33-e2) 

^{C^,[{t^Pa ■ e3)inpg ■ e^){xggf}a{T) - Xma{T)) 

+{lJ'fa ■ e3)(M/a • e4)XaaPa{T)] 

+C3[((M/^ • e3)(M/a • 64) - (Mas • e3)(M/35 " e4))x^a^a(T)]}, (11) 
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• e3)(M//3 • e4)x^/3^/3(T)] 
-Cf,Q(M/39-ei)(Ma9-e2) 

(T)) 

+Cf3[((M/a • e3)(M//3 • 64) - (M/3s • e3)(Mag " e4))Xa/3a/3(r)]}, (12) 



S{UJpg,T,UJpg) = -Cf^Cf,(M^g • ei)(/i^g • 63) 

x{Cf3[(A^/3<y ■ e3)(A^/35 ■ ei){xggi3i3{T) - 1 - XlsmiT)) 

+ {lJ'fa- e3)(M/a • e4)Xaa^^(7')] 

+C3[((M//3 • e3)(M/a • 64) - (Mag ' e3)(M^g " e4))xpaMT)]} 

x{C<^3[(^/39 • e3)(/i/39 • e^){Xggal3{T) - Xpfia/^iT)) 
+ {Hfa ■ e3)(At/a • e4)Xaaap{T)] 

+C^3[((/^//3 ■ e3)(M/a ■ 64) - {flag ' 63) (^^g " e4))X/3aa/3(T)] }. (13) 

Here, the expressions have been obtained using the rotating-wave approximation, as well as the 
assumption of no overlap between pulses, fipq — is the transition dipole moment between 
states p,q e {g,a, f}. We have rescaled the spectra amplitudes to eUminate the details of the 
Uneshape by multiplying them by the dephasing rates of the optical coherences in the coherence 
and echo intervals, 

S{U!pg,T,U!qg) = TgpTqgS{U!pg,T,(jJqg). (14) 

The coefficient C^. is the amplitude of the i-th pulse at the frequency Upg, 

18 



Echo 
time 

Waiting 
time 



71 



Coherence 
time T 



CO 



CD 



CO 



ag 



a(5 

aa 



ga 



gg 



CO 



ag 



CO 



ag 



CO 



ag 



ag 



a(3 

aa 



ga 



gg 



CO 



Pg 



CO 



ag 



FIG. 5. Dual Feynman diagrams that account for the population to coherence transfer terms Xai3aa(T) in 
quantum process tomography. 



CP = --v/2^e-"'('^™-^')'/2 (15) 

with A being the strength of the pulse and a the width of the Gaussian pulse in time domain. Also, 
Cj is the polarization of the i-th pulse. Both C^, and are experimentally tunable parameters for 
the pulses. 

Whereas Equations (14) and (15) presented in pT| correspond to a single value of r and t. 
Equations ( [TO] ), ( [TT] ), ( 12 1, and ( 13 1 stem from Fourier transform of data collected at many r and t 
times. Therefore, in principle, a 2D spectrum provides a more robust source of information from 
which to invert xiT) than in the suggested ID experiment. The displayed equations, albeit lengthy, 
are easy to interpret. For instance, consider the term which is proportional to Xapaa{T) in Equation 



(10), which stems from the Feynman diagram depicted in Fig. 5. As expected, it consists of a 
waiting time where the initially prepared population |a;) (q;| is transferred to the coherence ja) 
This waiting time is escorted by a coherence \g){oL\ oscillating as e*^~*'^9°~^9")'^ which evolves 
during the coherence time and another set of coherences and \(y){g\ which evolve during 

the echo time as e^"*'^^'^"'"/'^)* = e^^^'^^f^^^s)*. These two intervals correspond to the diagonal peak 
located at {uJag, (^ag)- Other processes that exhibit oscillations at those two respective frequencies 
appear as additional terms in the equation corresponding to that particular peak. 
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FIG. 6. Transfer of population in eigenstate |a)(a| to other populations and coherences in the eigenbasis 
of the single exciton Hamiltonian. The hierarchy equation of motion approach is used for a dimer system 
based on the parameters of the site 1 and site 2 subsystem of the Fenna-Matthews-Olson complex. Popu- 
lation in |a)(a| decreases iXaaaa{T), purple) and is transferred to |/3)(/3| iXi3i3aa{T), blue). Emergence 
of coherence from the initial population occurs in this model 0i-{Xa/3aa{T)} , yellow and Q{xaaaa{T)}, 
green). 

In Ref. dm, we showed that there are sixteen real valued parameters of which need to 
be determined at every value of T in order to carry out QPT of the single exciton manifold of a 
heterodimer. For an illustration, we shall describe how to obtain the elements XijaaiT). These 
quantities are shown in Fig. 6 and have been computed using the Ishizaki-Fleming model, with a 
bath correlation of 150 fs [fTOll . They display rich and nontrivial phonon-induced behavior, such 
as the spontaneous generation of coherence from a population in an eigenstate of the excitonic 
Hamiltonian, and therefore, is a very good example of how QPT provides access to this nontrivial 
information via the repetition of a series of 2D PE experiments. For this particular set of x(T) 
elements, we shall exploit the waveform of the pulses but not their polarizations, and for sim- 
plicitly we will assume the polarization configuration xxxx for each of the pulses including the 
heterodyning. 



Consider the possibility of using pulses with carrier frequencies centered about Uag and 
respectively, and such that their bandwidth is narrow enough that the pulse centered about uiag 
has negligible component at cu^g and vice versa. Then, we can carry out an experiment such that 
iTiFT; JT^i^ \r^\ ^ (experiment 1) for all i and notice that the diagonal peak at {uiag,u}ag) 

\C^-^^\ {C^^l I'-^t^sl 

reduces to: 
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(S'(a;«,T,a;„))^^^^ = -^^.^^^^^((Mas ■ ei)(Ma3 ■ e2)[(M/a ■ e3){nffs ■ 64) (16) 

(T) 

which implies that its evolution with respect to T directly monitors the transfer of the popula- 
tion prepared at |«)(«| to the coherence at Here, {■)xxxx denotes an isotropic average of 
the experiments performed with the xxxx polarization configuration. Xapaa{T) can be directly 
obtained if information of the dipole moments is known in advance. As can be checked easily, 
Xai3aa{T) = {xi3aaa{T))* Can, in principle, be also obtained directly from an experiment where 
^ 1 for all i (experiment 2) and monitoring {S{uJa,T,ujp))xxxx- Redundant measurements 
can be used as ways of effectively constraining the QPT. 

Similarly, the transfer from to other populations can be extracted by monitoring 

{S{ijJa,T,oja))xxxx in experiment 2 and {S{oJa,T,ujp))xxxx in experiment 1. These two linearly 
independent conditions are enough to extract Xggaa{T), Xaaaa{T), and XpfSaaiT), since there is 
a third independent condition based on trace preservation which reads Xggaa{T) + Xaaaa{T) + 

XliPaa{T) = 1. 

It is now important to verify whether the suggested experiments are feasible. In order to ensure 

\C~^ I 

conditions of the form Hft ^ 1, we need a ~ — - — ~ 470 fs, that is, the pulse needs to be 
long enough to guarantee the narrow band condition. This requirement, although attractive from a 
pedagogical standpoint since it yields block diagonal sets of linear equations, is a nuisance from a 
practical perspective, as decoherence mechanisms might be in the same timescale and might not be 
'seen' with such long pulses. However, the only essential requirement is a toolbox of two different 
waveforms for the pulses. A more sensible choice is a set of pulses centered about cOag and 
respectively, but having a ~ 100 fs. By carrying out 8 experiments alternating the two waveforms 
in each of the three pulses, each of the terms in Equations (10), ( [TT] ), (12), and (13) which are 



proportional to Cl_^Cl^C^^ for i,j,k G {a, j3} may be inverted to yield the block diagonal set of 
equations discussed above. 

In summary, we have presented three different tools for unraveling the role of quantum coher- 
ence in biological systems: a) techniques for obtaining the contribution of quantum coherences to 
biological processes; b) a microscopic simulation approach to explore the dynamics of these sys- 
tems by direct simulation; and finally c) a new theoretical proposal for an experimental procedure 
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that provides detailed information about the quantum procesess associated with energy transfer in 
the ultrafast regime. We believe that ultimately, a combination of these three techniques and tools 
from other groups will be collectively required to make definitive conclusions about the role of 
quantum coherence in photosynthetic complexes. 

We are thankful for financial support from the Department of Energy Frontier Research Center 
for Excitonics, the Army Research Office and the Sloan and Camille and Henry Dreyfus founda- 
tions. 
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